% This program generates artificial output data, derives the PVM-consistent
% current account and then:
% - estimates the model predicted CA series and compares it to the
% model-consistent simulated series in terms of correlation coefficients
% and variance ratios
% - tests the model both in terms of the regular Wald test of the cross-equation
% restrictions, the Wald test "linearized", and the F-test.


clear all;

lags=4; %number of lags in the VAR
r=0.01; %assumed real interest rate
initial=500; %number of observations to be weeded out to exclude effect of initial conditions
t=120; %sample size
T=initial+t;
nloop=10000; %number of data samples

volra=[]; corr=[]; test=[]; stat=[]; statexact=[]; statlin=[]; testlin=[]; testexact=[]; fcrit=[]; waldcrit=[]; kcoef=[];


% we will now generate nloop samples of T observations for Y using the estimated VAR for each country,
% and btain T observations of model-consistent CA.
% After discarding initial observations, we will use DY and the model-consistent CA to obtain the model predicted series

var = [-0.868017921	-0.804691022	-0.878794632	0.086874109	0.163174125	0.025932876	0.327579693	-0.502387565
      -0.106112655	-0.180637003	-0.20025133	-0.078293612	0.637950532	0.146085012	0.129228775	0.044684862]; %this is the matrix of VAR coefficients estimated from the country
companionvar = [var(1,:)
                1 0 0 0 0 0 0 0
                0 1 0 0 0 0 0 0
                0 0 1 0 0 0 0 0
                var(2,:)
                0 0 0 0 1 0 0 0
                0 0 0 0 0 1 0 0 
                0 0 0 0 0 0 1 0]; % companion matrix of var
            
loopcount = 1; 
for loop = 1:nloop;

    
dny=[];x=[]; u1=[]; u2=[]; ca=[]; capred=[];
z2=[]; lhs2=[]; x2=[]; x2p=[]; phi2=[]; res2=[]; x1=[]; phi1=[]; res1=[]; lhs1=[];
V11 =[]; V12=[]; V13=[]; V14=[]; V=[]; k=[];  

% we will now create artificial data using the data generating process

dny(1:lags)=1; x(1:lags)=1; ca(1:lags)=0;

sigma=[0.018712003 0.009110168
0.009110168	0.011825209];
u=mvnrnd([0 0],10*sigma,T);

for i=lags+1:T; 
    
    dny(i)=var(1,:)*[dny(i-1) dny(i-2) dny(i-3) dny(i-4) x(i-1) x(i-2) x(i-3) x(i-4)]'+u(i,1); %this generates data for dny according to the estimated VARs
    x(i)=var(2,:)*[dny(i-1) dny(i-2) dny(i-3) dny(i-4) x(i-1) x(i-2) x(i-3) x(i-4)]'+u(i,2);
    ca(i) = -[1 0 0 0 0 0 0 0]*(1/(1+r))*companionvar*inv(eye(2*lags)-(1/(1+r))*companionvar)*[dny(i) dny(i-1) dny(i-2) dny(i-3) x(i) x(i-1) x(i-2) x(i-3)]'; % %this is the model consistent current account
    

end;

z2=[dny(initial+1:T)'-mean(dny(initial+1:T)') ca(initial+1:T)'-mean(ca(initial+1:T)')];

% we are now going to estimate two systems
% 1) regress ca(t)-(1+r)ca(t-1)-dny(t) on the information set at time t-1,
% which gives us the 'exact' test
% 2) a VAR on ca and dny to obtain the predicted series and do the Wald
% test of the cross-equation restrictions of the model

for j=1:lags; x2=[x2 z2(lags+1-j:t-j,1)]; x2p=[x2p z2(lags+1-j:t-j,2)]; end; x2=[x2 x2p]; 

lhs1 = z2(lags+1:t,2)-(1+r)*z2(lags:t-1,2)-z2(lags+1:t,1); 
x1=x2(1:t-lags,:);                                                                                                                                                                                                                                        
phi1 = lhs1'*x1*inv(x1'*x1);                                                     
res1 = lhs1-(x1*(phi1'));
sig1 = (res1'*res1)/(size(res1,1)-2*lags);

lhs2 = z2(lags+1:t,:);  
phi2 = lhs2'*x2*inv(x2'*x2);
res2 = lhs2-(x2*(phi2'));
companionphi2 = [phi2(1,:)
                1 0 0 0 0 0 0 0 
                0 1 0 0 0 0 0 0
                0 0 1 0 0 0 0 0
                phi2(2,:)
                0 0 0 0 1 0 0 0 
                0 0 0 0 0 1 0 0
                0 0 0 0 0 0 1 0]; % companion matrix of phi2

% This will measure the highest eigenvalue of the VAR matrix
eigen=eig(companionphi2);
for i=1:8;
if real(eigen(i))==eigen(i);
eigen(i)=eigen(i);
else;
eigen(i)=NaN;
end;
end;
eigmax(loopcount,1)=max(eigen);  
eigmax2=eigmax(~isnan(eigmax));

% we are now constructing V, the variance-covariance matrix of VAR
% coefficients

V11 = (res2(:,1)'*res2(:,1))/(size(res2,1)-2*lags)*inv(x2'*x2);
V12 = (res2(:,1)'*res2(:,2))/(size(res2,1)-2*lags)*inv(x2'*x2);
V21 = (res2(:,2)'*res2(:,1))/(size(res2,1)-2*lags)*inv(x2'*x2);
V22 = (res2(:,2)'*res2(:,2))/(size(res2,1)-2*lags)*inv(x2'*x2);
V = [V11 V12;V21 V22]; diag(V);

k = -[1 0 0 0 0 0 0 0]*(1/(1+r))*companionphi2*inv(eye(2*lags)-(1/(1+r))*companionphi2);
kcoef(loopcount,:)=k;

capred(1:lags) = 0;

for i=lags+1:t;
    
    capred(i) = k*[z2(i,1) z2(i-1,1) z2(i-2,1) z2(i-3,1) z2(i,2) z2(i-1,2) z2(i-2,2) z2(i-3,2)]';
    
end;

% this constructs the Jacobian of derivatives of K with respect to the VAR
% parameters

for i = 1:4*lags;
        a=[]; varvec=[]; kp=[];
        a = reshape(transpose(phi2),1,4*lags);
        a(i) = a(i)+0.000001;
        varvec = [a(1:2*lags)
                  1 0 0 0 0 0 0 0
                  0 1 0 0 0 0 0 0 
                  0 0 1 0 0 0 0 0
                  a(2*lags+1:4*lags)
                  0 0 0 0 1 0 0 0
                  0 0 0 0 0 1 0 0
                  0 0 0 0 0 0 1 0];
              
        kp = -[1 0 0 0 0 0 0 0]*(1/(1+r))*varvec*inv(eye(2*lags)-(1/(1+r))*varvec);
        J(1:2*lags,i) = (transpose(kp)-transpose(k))/0.000001;

end;

% We will now construct the 'exact' test...

statexact(loopcount,1) = ((lhs1'*lhs1-res1'*res1)/(2*lags))/((res1'*res1)/(size(res1,1)-2*lags)); %this is the F-statistic
fcrit(loopcount,1)=fcdf(statexact(loopcount,1),2*lags,size(res1,1)-2*lags);

% and the wald test
stat(loopcount,1)=(k-[0 0 0 0 1 0 0 0])*inv((J*V*J'))*(k-[0 0 0 0 1 0 0 0])';
waldcrit(loopcount,1)=chi2cdf(stat(loopcount,1),2*lags);

% and the linearized wald test

R=[-eye(2*lags) eye(2*lags)];

statlin(loopcount,1)=(R*[phi2(1,:)';phi2(2,:)']-[0 0 0 0 1+r 0 0 0]')'*inv(R*V*R')*(R*[phi2(1,:)';phi2(2,:)']-[0 0 0 0 1+r 0 0 0]');
waldlin(loopcount,1)=chi2cdf(statlin(loopcount,1),2*lags);

if waldcrit(loopcount,1) <= 0.95;
    test(loopcount,1) = 1;
else;
    test(loopcount,1) = 0;
end;

if waldlin(loopcount,1) <= 0.95;
    testlin(loopcount,1) = 1;
else;
    testlin(loopcount,1) = 0;
end;

if fcrit(loopcount,1) <= 0.95;
    testexact(loopcount,1) = 1;
else;
    testexact(loopcount,1) = 0;
end;


% this computes the variance ratio and correlation coefficient between predicted and actual

z=corrcoef(capred(lags+1:t),z2(lags+1:t,2)); 
corr(loopcount,1)=z(1,2);
volra(loopcount,1)=(std(capred(lags+1:t))/std(z2(lags+1:t,2)))^2;

loopcount=loopcount+1;
end; % end of loops
